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PREFACE 


This final report summarizes the. work completed under Contract Ho. 
NASi-12726 towards the development of models and associated computer codes 
for the analysis of three-dimensional supersonic nozzle-exhaust flowfields. 
The contract was monitored by Mr. Manuel Salas of NASA Langley who pro- 
vided comparison flowfield 1 results and constructive suggestions during the 
course of this effort. The authors additionally acknowledge the benefit 
of many fruitful discussions with the late Dr. Antonio Ferri and the pro- 
gramming assistance provided by Mr. Paul Kalben in the development of com- 
puter codes BIGMAC and CHAR3D. 
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I . INTRODUCTION 


This report summarizes work accomplished under Contract No. NAS1-12726 
towards the development of computational procedures and associated numerical 
codes for three-dimensional nozzle-exhaust flow fields. The flow fields con- 
sidered are those associated with airbreathing hypersonic aircraft which re- 
quire a high degree of engine/airframe integration in order to achieve opti- 
mized performance. The exhaust flow, due to physical area limitations, is 
generally underexpanded at the nozzle exit; the vehicle afterbody under- 
surface is used to provide additional expansion to obtain maximum propulsive 
efficiency. This results in a three-dimensional nozzle flow, initial ized at 
the combustor exit, whose boundaries are internally defined by the under- 
surface, cowling and walls separating individual modules, and externally, by 
the undersurface and slipstream separating the exhaust flow and external 
stream. A typical exhaust nozzle is depicted in Figure (t), characterized 
by multiple rectangular nozzle modules. 

The numerical models, developed in this analysis address the following 
characteristic features of these exhaust flows: 

(1) The flow properties at the combustor exit are highly nonuniform. 
Burning and mixing in the combustor yield regions of highly varying 
composition, temperature, and stagnation properties, in addition, 
shock waves are produced in the vicinity of the injectors. Although 
the strength of these waves decays rapidly as they propagate through 
the burner, they are generally present at the burner exit and must 
be accounted for. 

(2) The exhaust gas mixture consists of hydrogen-air combustion products 
and significiant burning may still occur in the initial regions of 
nozzle expansion. 

(3) The flow field geometry is quite complex. The engine modules con- 
sist of multiple u.rfaces with sharp interior corners, and flow 
fences, to contain the external exhaust flow, may be present. 
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(4) The interior nozzle flow field is dominated by complex wave inter- 
actions with waves generated and reflected off multiple surfaces. 

!n addition, sharp interior corner regions must be accounted for. 

(5) The nozzle exhaust flow interacts with the nonuniform vehicle ex- 
ternal flow field. This complex interaction for underexpanded ex- 
haust flows results in an expansion system propagating toward the 
vehicle undersurface from the cowl trailing edge and a spanwise 
expansion generated by the sidewall interaction. An underexpansion 
shock propagates outward into the nonuniform vehicle external flow, 
and the exhaust and external flow are separated by a plume boundary. 
In addition, pressure and flow deflection mismatch between adjacent 
modules may occur, resulting in a spanwise multiple shock system. 

To accommodate the varied, complex geometric configurations entailed in 
this analysis, a reference plane approach has been utilized, with respect to 
several coordinate systems. This approach involves the definition of a re- 
ference plane system in which the three-dimensional volume under consideration 
is spanned by an appropriately selected series of planes which intersect the 
bomdaries of the considered volume. The equations of motion within the re- 
ference planes are expressed in a quas i-streaml ine coordinate system, where 
quasi-streamlines are the projections of the actual stream surfaces onto these 
reference planes. Such a system accommodates the calculation of highly rota- 
tional, variable composition flow fields by minimizing streamline interpola- 
tion procedures which can produce significant errors, as discussed by Sedney 
in Reference (1) . 

Reference plane systems in cartesian, line source and cylindrical coordi- 
nates are illustrated in Figures (2), (3) and (4) respectively. The configura- 
tion of the reference planes is chosen to best accommodate the overall flow 
field geometry by having primary flow variations occur within the -eference 
planes. A more complex flow field situation is depicted in Figure (5) for 
the flow downstream of the cowl exit. For this calculation, a combination of 
several reference plane systems is employed and provisions are included in 
the numerical codes for automatic switching from one system to another as the 
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. character of the boundary surfaces changes. The reference plane system also 
caters to the usage of reference plane characteristics at all boundary points. 
This approach is generally recognized as the most accurate boundary calcula- 
tional procedure (Reference 2). However, it proves cumbersome when employed 
In conjunction with nonreference plane networks due to the complex inter- 
polation procedures then required. 

The reference plane characteristic technique has been widely used for 
the calculation of three-dimensional supersonic flow fields, and the authors 
had previously developed a program employing this technique for the calcula- 
tion of nozzle exhaust flow fields (References 3 and 4), which is in current 
usage at NASA Langley Research Center (References 5, 6 and 7). That program, 
as well as most reference plane characteristic (refchar) codes in common 
usage (References 8 and 9), employs an inverse scheme wherein interpolations 
are performed to obtain data at the intersection of the quasi -characteristics 
with a noncharacteristic initial data surface. Comparisons of such refchar 
codes with shock capturing finite difference codes (Reference 10) have led to 
the general cone 1 us ion. that difference codes are better able to analyze com- 
plex flow fields with multiple secondary shocks. From experience gained with 
the authors’ original refchar code, it was felt that the inability to success- 
fully analyze such flow fields was primarily due to the inverse interpolation 
procedures employed. Such procedures tend to ignore the presence of weak waves 
by allowing the quasi-characteristic lines to arbitrarily cross each other. 

The numerical diffusion associated with these interpolations can become sig- 
nificant, particularly when the local Courant number (ratio of overall marching 
step to local maximum allowable marching step) is much less than one. The 
smearing of these weak waves is enhanced by resorting to higher order interpola- 
tions on the initial data line. 

To treat complex multiwave flow fields and still retain the advantages that 
reference plane methods afford, two new numerical codes have been developed. 
Program CHAR3D is based upon a quasi-characteristic calculational procedure 
employing a novel wave preserving network, as compared to previous standard in- 
verse networks. In addition, a nonisentropic pressure-density relation applied 
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along streamlines, permits the estimation of shock entropy losses while, 
the usage of conservation variables in the construction of cross-flow deriv- 
atives has facilitated the analysts of flow fields containing cross stream 
discontinuities. Program BIGMAC is a reference plane finite difference code 
utilizing the described quasi-streaml ine grid within reference planes. Shock 
eapturing capabilities are provided via the usage of conservation variables 
in conjunction with a one-sided difference algorithm. 

While the numerical algorithms and logical procedures employed in CHAR3D 
and BIGMAC differ, both codes employ a newly developed geometry package (Re- 
ference 11) for a description of boundary contours, calculate boundary points 
by reference plane characteristic procedures, and incorporate the same thermo' 
dynamic fits far describing the hydrogen-air gas mixture in chemical equilib- 
rium. 


The governing flow field equations are presented in Section M while 
computational procedures for CHAR3P 'nd BIGMAC are presented in Section III. 
Section IV contains a description of rumple calculations performed with both 
these codes. Conclusions drawn from this study and recommended procedures 
in extending this analysis are presented in Section V. For completeness, a 
sunanary of the equilibrium curve fits are presented in Appendix I while a 
description of the geometry package is presented in Appendix it. 


2®^ page & 

DU POOR QUALITY 


9 


II. GOVERNING EQUATIONS 


Ac Characteristic Analysis (CHAR3D) - The equations governing the steady 
flow of an inviscid gas mixture in chemical equilibrium may be written: 


Continuity 

V» (pV) =0 (1) 

Momentum 

p(V«V) V + VP = 0 (2) 

Energy 

V*VH * 0 (3) 

Equivalence Ratio Constancy Along Streamlines 

V*V$ « 0 (4) 

This system is supplemented by the relation 

V-VS =0 (5) 

expressing constancy of entropy along streamlines in continuous regions of the 
flow field. The equation of state may be written 


(IE) 

V, 



2 

a 


where curve fits for the isentropic exponent F, from Reference (4), 


( 6 ) 


r - f (h, p, $) 


(7) 


are described in Appendix I. The continuity equation, employing Equations (5) 
and (6) may be written: 

V*VP + a 2 pV*V * 0 (8) 


The scalar forms of the modified continuity equation (Equation 8) and the 
momentum equations (Equation 2) in general orthogonal coordinates are: 
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Modified Continuity 


u w • ' tf F w P * 

2 x t X 1 

pa (jj — + ij-*) + -f; + — r~fr~jr*- th.h (-^* P + pv ) + 

n t n 3 n 1 n 3 n 1 h 2 h 3 1 3 a 1 x 2 x 2 


h 2 ■ h_ pu + h.hg pw] 
X 1 x 3 


x^ Momentum 


“ U t « u + 1 p = 

ph^ 


h 2 U>c 2 + h i h 2 


1 


Xg Momentum 


+i 2- v + _1_ p = - ^ v 

h 1 X 1 h 3 x 3 0h 2 x 2 h 2 x 2 h l h 2 2 x. 


vw . 
h 2 h 3 2 x. 


Momentum 


JL„_ +«„ + 4 . P =. » w + * h , 

h l "1 h 3 x 3 ph 3 x 3 h 2 x 2 h 2 h 3 2 x- 


The Identification of metric coefficients and coordinate directions for the 
three coordinate systems under consideration is provided below. 


System 

X 1 

X 2 

X 3 

h i 

h 2 

h 

h 3 

Cartesian 

X 

y 

z 

i 

l 

1 

Cylindrical 

X 

e 

r 

i 

r 

1 

Line Source 

r 

e 

z 

i 

r 

1 


Note that h 2 is unity for the line source system and zero for other systems 




while hj is unity for the cylindrical system and zero for other systems. 


x, 
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These terms will be replaced by the Indices and where * 1 for the • 
line source system and J 2 » 1 for the cylindrical system. Both indices are 
zero otherwise. 

In writing Equations (9) ~ (12), the source terms and terms Involving 
derivatives normal to the* reference planes - constant have been put on 
the right-hand side. Then, the left-hand side of these equations corresponds 
to that of the two-dimensional system in the x^, x^ plane. Reference plane 
characteristic relations are readily obtained by algebraic manipulation of 
the modified continuity equation and the x^ and x^ momentum equations. 

These relations are listed below, where the velocity vector 

A A A 

V*&u{ + v I + w i 

X 1 x 2 x 3 

is expressed in terms of its magnitude in the reference plane, q, quasi- 
st.eamline angle, and cross flow angle, employing the geometric rela- 
tions 

o •) ^ 
q * (u 4 - w ) 

$ “ tan 1 (w/u) 

^ * tan" 1 (v/q) 

as depicted in Figures (6), (7) and (8). Along the reference plane char- 
acteristics 

+ ^ ^3, = M 2 costj>s T nd»±S 

dx 1 M 2 eos 2 (<f>-1) 

the compatibility relation may be written: 

d$ ± i d(£n P) = F" dx 

rtr 


(13) 


(14a) 

(14b) 

(14c) 


(15) 

( 16 ) 
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V a U i x + V 
Q = V cos i|> 

u = Q cos <p 

v - Q tan ij> 

w = Q sin 4> 


- 2 ^ 

Q = (u 2 +w 2 ) 

$ * tan"*(w/u) 

if» = tan” 1 (v/Q) 


FIGURE 6. VELOCITY VECTOR IN CARTESIAN SYSTEM 



V 


FIGURE 7 



Q = 
u = 
v = 
w = 


, i x + v i 0 
V cos 4) 

Q cos $ 

Q tan 4 > 

Q sin 4> 


+ w i r 

Q = (uW)** 

4 > = tan" 1 (w/u) 
4 j = tan" 1 (v/Q) 


VELOCITY VECTOR IN CYLINDRICAL SYSTEM 



Q “ V cos 0 - (u^+w^J 

u = Q cos 4> 4> = tan“*(w/u) 

v = Q tan ij> i|) = tan“*(v/Q) 

w = Q sin <J> 


FIGURE 8. VELOCITY VECTOR IN LINE SOURCE SYSTEM 



where 


H 2 « q 2 /a 2 , e 2 - M ? - 1 


and 


F + » (sin# - A + cos<fi) [(tantfi) + Un P) ] 

x 2 -2 
f 4. 4 . 2 

- <ji taniji {ccs$ + A - sin 4 >) + (J 7 - A“J 9 ) tan 

Xg * * 


(17) 


dx 


1 


dx = — j- + d in Xj 

*8 2 

The streamline projections onto the reference planes 

dx„ 

A, 


‘SL “ d^7 * tan * 


( 18 ) 


are also characteristic directions for this system. Along the quasi-stream 
Ag^, the relation between cross flow angle, ip, and pressure, P, may be written 


d(tamb) = d(&n P) + G dx (19) 

rM^ 

where Un P) x 

G = -J— - t 5 — - + tan*(tan*) 

COS9 ^2 

(20) 

4* tanif)(l + tan 2 ^) (4^cos«ji + J^s i n«f») 1 


B. Finite Difference Analysis (BIGMAC) - The analogous inviscid flow 
equations in conservation form may be written: 

Continuity 


V-(pV) - 0 

Momentum 

V.(pW) + 7P = 0 (21) 
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Energy 


7-{pVH) * 0 


( 22 ) 


Equivalence Ratio Constancy Along Streamlines 


7*{pV$) =* 0 ' (23) 

t * “ . 

The scalar form of these equations in the x^ , x^, x^ reference plane system 
may be written: 


+ F + G + H 


where for k * 1 to 6 


(24) 





In performing the numerical integration of this system of equations, a quasi- 
streamline grid network is employed wherein one follows the projections of 
streamlines on the reference planes and one to one correspondence is made be- 
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tween "corresponding" streamlines on adjacent reference planes. Then Equa- 
tion (24) may be written 


+ F a, + G 

X 1 X 2 


y + H - r~ - — E 
x^ hj u Xj 


"2 - 
tana t— F 

h 3 *3 


a. 


where 3/3x^ denotes the partial derivative along the quasi-streaml ine, in the 
reference plane « constant while 3/3x 2 denotes the partial derivative along 
the line connecting "corresponding" streamline points on adjacent reference 
planes at the marching station, = constant. Thus, 


(_L_) 


= W 

X 2’ X 3 1 x 2 ,n 


if- 

~iH (JL) 

h^ u '3Xj 




and 


«*£ 


x r x 3 


(-*-) 

3x„ 

* . z 


hf “» <%> 


x 1 


x 1 ,x 3 


where tana is the slope made by the line connecting "corresponding" grid points 
in adjacent reference planes with respect to the x^ coordinate direction, at 
x.j « constant. 
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III. COMPUTATIONAL PROCEDURES 


A. Reference Plane Grid Networks - Both CHAR3D and BIGMAC employ refer- 
ence plane grid networks following the trace of streamline projections onto 
the reference planes as illustrated in Figure (9). Note that the grid network 
utilized for CHAR3D (Figure 9a) differs from that of the author's previous 
characteristic method (References 3 and k) as well as from that of other 
popular reference plane characteristic methods (References 9 and 10). Pre- 
vious methods employ a standard inverse scheme (i.e.. Figure 9b) wherein 
Interpolations are performed along a non-characteristic Initial data surface. 
Such Interpolation procedures ignore the advantages affordable by a reference 
plane characteristic method by failing to utilize information contained along 
already calculated quas i -characteristic surfaces. By interpolating along such 
surfaces, rather than along non-characteristic initial data surfaces, two 
distinct advantages accrue: 

(1) The wave nature of the flow field within tne reference planes is 
preserved. • . 

(2) Linear interpolation procedures along such surfaces is compatible 
with a scheme of second order accuracy (see Reference 12, appendix). 

The application of linear interpolation procedures along a non-characteristic 
initial data surface does not produce results accurate to second order, while, 
implementation of higher order interpolative procedures tends to result in ex- 
cessive numerical diffusion (see Reference 1). 

While both CHAR3D and BIGMAC empioy the same quasi-streaml ine reference 
plane network, the computational sequence differs substantially. Let I desig- 
nate a grid point in reference plane J; for the calculation of an internal 
module such as that in Figure (2) each reference plane J contains IMAX grid 
points where 1 = 1 represents the lower boundary surface (i.e., vehicle 
undersurface) while i = IMAX represents the upper boundary (i.e., cowl), in 
this example J = 1 is a plane of symmetry while J * JW indicates the module 
sidewall. In CHAR3D, the lower boundary point i = 1 is first calculated for 


-a 



□ WAVE GRID POINT 
O STREAMLINE GRID POINT 


(a) Reference plane grid network for CHAR 3D. 



0t>) Standard reference plane 
characteristic network. 



(c) Finite difference network. 


Figure 9- Reference plane networks. 
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a] I reference planes J = 1 to JW, then the grid points I * 2 for J * 1 to JW, 
etc., proceeding to the upper, boundary point I = IMAX. This provides a quasi- 
characteristic initial data surface comprised of the dawnrunning reference 
plane characteristics passing through the points I for J =* 1 to JW for the 
calculation of the points I + 1 for J = 1 to JW. 

Program BIGMAC proceeds in the opposite fashion. Predictor values are 
calculated for ail grid points I = 1 to IMAX in a particular reference plane 
J starting with J = 1 and proceeding to J = JW. The process is repeated for 
corrector values. 

Both programs employ internal disc storage to provide for ti.e usage of a 
large number of grid points without exceeding the small core memory alloca- 
tions of the CDC 7600. Thus, CHAR3D provides storage for all reference 
plane locations J = 1 to JW (JW <_ ^3) at 5 levels of I while BIGMAC provides 
storage for all reference plane grid points I = 1 to IMAX (IMAX <_ 40) for 10 
reference planes J. In both programs, reference planes are deleted or added 
in the proximity of sidewalls according to the following criterion. Let A 

represent the spacing between reference planes and A the spacing between the 

w 

last reference plane and the sidewall. Then, the reference plane adjacent to 
the wall is deleted when A w < .6 while a plane is added between the last re- 
ference plane and sidewall when A > 1.6. 

w 

B. Interior Point Calculational Procedure - Properties are desired at 
the grid point (I,J,K) shown in Figure (10) for a cartesian system. The 
allowable step size Ax is determined by satisfying the CFL condition. For 
BIGMAC, this requires that the intersection of the Mach cone from (I,J,K) 
with the initial data surface fs?l 1 s within the numerical domain as depicted 
(i.e., the quadrilateral (I, J + 1), (I - 1, J) , (l, J - 1), (l + 1, J)). 

Note that the effective numerical domain for the characteristic calculation 
includes, the points I + 1 and I - 1 on planes J - 1 and J + 1; hence, a 
larger step may be taken with CHAR3D (Ax fHAR3D y 2 a * B!GHAC )- 
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CHAR 3 D 


Referring to Figure (11), for a cartesian system, the grid point T is 
located along the quasi -streaml ine through grid point 1 by the difference 

approximation to Equation (18). 

» 

Z T J = Z ! j + tan< ! , i j + tantj>j j) Ax (26) 

where a =* 1, b - 0 in the predictor step and a = yy, b * y in the corrector 
step. In this new wave preserving network, the calculation proceeds upward 
from the lower boundary where points (T - 1 , J) are calculated for al 1 re- 
ference planes J to second order prior to calculating points (T,J). In addi- 
tion to the standard initial data array (the points (l,J)), an extra array 

a — 

(|,J) is required. To calculate properties at (I,J)» the standard initial 
data grid in the reference plane (I - 1, J) , (1,J), and (I + 1, J) is employed 
to calculate the forcing function terms involving derivatives normal to the 
reference planes. Properties are known at points H ^ , G^, and 1-1 from the 
calculation of point (T.- 1, J) to second order. 

4 * — 

Point A is located between and on the quasi -characteristic X“(Al) 
where X ± Is defined in Equation (15). All properties (including forcing func- 
tions) are obtained via linear interpolation between and G^ . Then, is 
located between ! and I + 1 such that the downrunning quasi -characteristic from 
G^ (or B) passes through (T,J). To first order, properties at (l,J) are cal- 
culated using points B and A, where Pg and <J> B are determined using compatibility 
relations (Equation 16) along IB and H^B. 


Then Py and tpj are calculated employing the compatibility relations 


and 


(#y " <1 > ft ) + 


(+ T - y - 


a (-%> + b (-2-) 

L rH A ™ T-l 


G 


Pt 


An (~) = (aF+ + bF*) Axy ft 
A 


(27) 


a (-^) + b (-^) 

™ B Tn y 


An (^-) =* (aF“+ bFy) Axy Q 
B 
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REFERENCE PLANE J (y ■CONSTANT) 



O STREAMLINE GRID ARRAY 

• INTERMEDIATE POINT 
(CHARACTERISTIC 
CALCULATION) 

□ EXTRA CHARACTERISTIC 
ARRAY (INTERPOLATED 
ALONG CHARACTERISTIC) 


Figure 11 . CHAR3D interior point grid. 
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Remaining properties are determined at I via the following streami irte rela- 
tions: 


(taroi>)j = (tan4»)y + 


a ( *•£) + b ( tanjL) 

™ | rM T J 


P T 

An (J-) 
P l 


+ (aGj + bGj) Ax | j 


( 28 ) 


h T ' H l ' 


a (tSSt H ) + b t-ESSi H ) 

cos$ y' j cos<f> y y 


' ^ 
Ax 


I i 


(29) 


“ 4>j - 


cos(j> y j 

and in continuous regions of the flow 


a { 


tan£ * ) + b ( tan£ # } 

cos<J> y y 


.'V 

Ax,! 


(p/p r ) Y = (p/p r ) I - 


a [“S- (P/P F ) J + b [|S0i (p/ p r )l 


COS<j> 


I 


cos^t 


I 


(30) 


Ax jT (30 


The following three parameter curve fits (based on data from Reference 13) 
are incorporated into this code and are described in detail in the appendix 
extracted from Reference (4). 


h » h (P,$,T) 
P “ P (P>$,T) 

r * r (p,$,t) 


(32a) 

(32b) 

(32c) 


The flow velocity is obtained via the relation 

Vy = -Ti [Hy - h (Py, $y, Ty)]* (33) 

where Ty is obtained via an inversion of Equation (32b) (with py, Py, and $y 

known) and h- is obtained employing Equation (32a). Then, Ty is obtained from 

2 

Equation (32c) and a - fyPy/py. 
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This calculation is performed for points 7 In al 1 reference planes to 
first order. Then, cross derivatives 3/3y are evaluated at 7 employing the 
relation 


/ 3f \ _ f 3f \ 

3y " (ay> " Lanct 3z 

7 x,z x,n x,y 


where 


(— ) 

v 3y' 


x,n 


Ay. Ay Ay 

„ J r_ / l _ 1 

l,J+1 y Ay 2 ' r 7,J ^Ay t Ay, 

(Ay 1 + Ay z ) 


i Ay_ 

r J " f 7,j-i { Zyf 


( 34 ) 


Ay t 



y 7,j-i 


Ay. 


VT 


I , J+1 


' y T,J 


tana - (f^) 

x,n 

f- - f- 

f) = JjJ Jzhl 
92 x,y Z 7,J " z 7 -1 ,J 

Derivatives are made the same way at the initial station !, except here 3f/3z 
is evaluated by 

(it) - ~ f i.j 

3z x,y " 2 -H,J - Z I,J 

CHAR3D» In addition to the centered difference algorithm described above, 
has the option of evaluating cross derivatives via an alternating one-sided 
difference algorithm. For this option, derivatives are evaluated as described 
in the section for BIGMAC. Cross derivatives are required for the variables 

p 

P, <J>, H, and P/p . In evaluating cross derivatives for P, <j>, and if>, 
conservation variables are employed as follows: 
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p y " F(3) y ‘ F(1) v y " v F(,) y 


■* 


w cosd) - u sin4> 

. y v I 

y q 

qv - (uu, + ww ) tarn); 




v . v y 


where 


E(3) - v E<1) 

V y eTH * 


E(2) - h„h,P - u El) 

- ZHv n 

e nr 


w 


E(ii) y - w E(l) y 

EOT 


( 35 ) 


( 36 ) 


The conservation variabl.es E(k) and F(k) are given by Equation [Zb). The 
use of conservation variables in the construction of these cross derivatives 
has tended to suppress oscillations that occurred when employing physical 
variables to difference across shock waves. However, the use of a one-sided 
difference algorithm in conjunction with ChAR3D has tended to produce spurious 
results in regions of large cross flow. 


In the characteristic reference plane algorithm, cross flow variations are 
expressed via the forcing function terms F” appearing in the right side of the 
compatibility relations (Equation 16). These terms are assumed to vary mildly 
within an integration step. When a one-sided algorithm is employed to evaluate 
cross derivatives in the vicinity of shocks, the values of the forcing function 
terms may vary greatly between the predictor and corrector steps. In addition, 
the numerical domain of dependence is somewhat vague for the characteristic 
reference plane approach in conjunction with one-sided differences, so that 
part of the problem encountered may be due to stability. The recommended ap- 
proach for evaluating cross derivatives in CHAR3D is to employ conservation 
variables in conjunction with a centered difference algorithm, although this 
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matter requires further study. 

In CHAR3D, secondary shocks are captured as rapid changes spread over ap- 
prpximately three grid points. The entropy change associated with these shocks 
is evaluated employing a nonisentropic pressure-density relation (illustrated 
here for a perfect gas). 

V*7 An (P/p Y ) o V • ~ ( 37 ) 

v 

For a shock of strength £ (pressure ratio across shock), this change is deter- 
mined employing the relation (for perfect gas) 


AS 

= An £ - y An 


r (y + 1 )g + (y - 1 ) 
l (Y “ + (Y + 1) 


] 


( 38 ) 


where AS is the entropy change along a streamline produced by the captured shock. 
This relation involves only the pressure distribution in the vicinity of the 
shock and is readily applied in regions of non interacting shocks as follows. Let 


F(e,r) 


(y + ik ± (y - i) 
(y ■ iTf + (y + i) 


Assume a shock is spread over the marching interval K = 1 to 6 (Figure 12) for 
a typical quasi-streamline. Then 1 represents free stream conditions for this 
shock. The entropy change in the interval K - 1 to K is then expressed by 


(f-) = (f) - <~) 

V K-1,K V I ,K V I , K-l 


= An 


’1 ,K 1 1 K- 1 


L C I,K-1 F 1,K 


where 


,K ” P K /P 1 


Then 


(p/p y ) k = 


m {e, °\ a Vi.k 


exp (~) 
U v 


( 39 ) 


K-l ,K 
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-Since the shock geometry does not appear in the entropy jump relation, the 
entropy rise associated with extremely complex three-dimensional shocks can 
be accurately obtained. Special provisions have been incorporated into the 
program for the computation of singular points at the juncture of intersecting 
shockwaves and/or shock reflection points. At such points the streamline 
undergoes a discontinuous pressure rise corresponding to that through both 
shock waves. If the shock intensities are different, an entropy discontinuity 
occurs separating the different zones, and a vortex of infinite intensity re- 
sults. Numerically, the entropy procedure described would predict an entropy 
rise associated with this pressure jump. Theoretically, this occurs in the 
limit of vanishing mass flow, while numerically the finite mass within this 
region would lead to unduly large entropy levels. Special coding has been 
incorporated at such singular points to suppress these "numerical" peaks. 


B1GMAC 


The HacCormack difference algorithm of Reference (1*») is employed for the 
calculation of interior grid points in BIGMAC. Referring to Figures (13)* 
(14) and (15) for grid index notation in the coordinate systems considered, 
this two step predictor-corrector scheme as applied to Equation (25) yields 
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= E, 


l,J U I,J 


Predictor Step 
h 
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The physical variables are obtained by the following iterative decoding 
procedure. A value of u is assumed. Then, 


3* 


(40b) 



p » E(l)/(h 2 h 3 u) 

P - (E(2) - E{1) u)/(h 2 h 3 ) 


(41a) j 

,• \ 

(4ib) 

. 4 


(hie) 

1 

1 

(4ld) j 
(4le) ] 

** ■ E(6)/E(l) (41 f) | 

h =» H - (u 2 + v 2 + w 2 ) (4lg) j 

•The value of h obtained in Equation (4lg) yields T via an inversion of Equa- ] 

tion (32a). Equation (32b) yields an alternate value of the density compared jj 

to that obtained in Equation (4la), The value of u is perturbed and the i 

procedure repeated until the two values of density agree to within a specified 
tolerance. A linear error extrapolation is employed to speed convergence. 


v * E(3)/£(1) 
w = E(4)/E(l) 
H » E(5)/E(1) 


Note that both codes additionally provide for the calculation of a uni- 
form composition, perfect gas mixture. For this option, the equilibrium 
sound speed, a, of Equation (6) is replaced by the frozen sound speed 


2 

a f 




m 


where the constant specific heat ratio of the frozen mixture y, replaces the 
equilibrium isentropic exponent T in all relations. The static enthalpy is 
then expressed by 


h (P,p) = P/p 


(43) 


and the iterative decoding procedure of Equations (4la-g) may be replaced by 
the direct determination of the u velocity component via solution of the 
quadratic (Reference 15) 


page to 
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m 


u - [- B + (B 2 - AAC) ]/2A 

where - - 

A => (y+1 ) /2y 

B » - E(2)/E(1) 

C * nl ( 2H - v 2 - w 2 ) 

C. Solid Surface Calculat ional Procedure - Surface geometry Is pre- 
scribed via discrete contour data for all continuous surfaces comprising the 
nozzle/afterbody configuration. Thus, for rectangular modules, as depicted 
in Figure (1), contour data is provided for four surfaces, namely, the two 
sidewalls, vehicle undersurface and cowl. These surfaces are fit via the 
method of Reference (11) based on the use of partial cubic splines; a sum- 
mary of this fitting technique is provided in Appendix II. Surface fitting 
Is done external to programs CHAR3D and BIGMAC via program FIT3D, (Reference 
It) which generates ordered surface coefficient arrays read in via tape to 
CHAR3D and BIGMAC and employed in these codes in conjunction with surface 
interpolation procedures’ also described in Appendix II. 

Both CHAR3D and BIGMAC employ reference plane characteristic procedures 
in the performance of a 1 1 boundary calculations. The calculational procedure 
for CHAR3D closely follows the step by step procedures detailed in Reference 
(4) . Referring to Figure (16), which depicts a lower boundary calculation in 
cartesian coordinates, CD is the intersection of the reference plane y = y^ 
with the lower surface z = f(x,y). In CHAR3D , the following iterative procedure 
is entailed: 

(a) 

(b) 

sfn<ji c = (f x ) cos«j» c + (f ) tamj» c (1»5) 


The cross flow angle iji is assumed equal to that at point D. 

c 

_ A 

The boundary condition V*n = 0 applied at point c yields the 
relation 
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(c) P c is obtained employing the compatibility relation (Equation 

27) along X using the above value of $ . 

c 

(d) An alternate value of ip c is obtained employing the normal mo- 
mentum relation (Equation 28) along CD. 

This new value of is then employed and steps (b) , (c) and (d) are re- 

0 

peated. This self convergent procedure generally requires only 2 or 3 itera- 
tions, although more may be required in regions of strong cross flow. Upon 
convergence, streamline relations yield remaining properties at c and the entire 
process is repeated for second order accuracy. 


In BIGMAC, this iterative procedure is eliminated by combining the normal 
momentum equation with the quasi-streamline momentum equation yielding the 
following system of equations for P , (w/u) and (v/u) in general coordinates 

C C C 

for the boundary x^ = f(x^> x 2 ) : 
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(f x > 
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and 


vP 


A 


6(1) ' G(1) 


1 h,h 3 (-/ + p» y ) + J, ^ - J 2 „ 2 


B - F(2) - u F(1) - J, 


v Fd) 


1 


C - F(4) y - w F(1) y - J 2 


E su - F(3) y - V F(1) y + J, 1121 + J 2 IIS! 


Note that in the above relations, v^ and P^ are evaluated in accordance with 
Equation (36). 


An important consideration in the boundary procedure for BIGMAC is the de- 
termination of the entropy at point C. The entropy change along the wall from 
D* to C is set equal to that along the streamline one mesh interval away from 
the wall. Thus, entropy changes associated with captured shocks as determined 
for interior grid points, are reflected in the wall point calculation. To 
accommodate this procedure, the lower wall point calculation must be deferred 
until after the calculation of the grid point one mesh interval from the wall 
and, the left sidewall calculations must be deferred until after the calculation 
of the adjacent reference plane. Then, predictor values are obtained at point C 

for P , (w/u) and (v/u) employing coefficients evaluated at the initial sta- 
CO c 

tion. Streamline relations yield H and $ and the density p f via Equation (39) 

o c c 

where the entropy change is evaluated as described above. The velocity magnitude 

V is obtained via Equation (33) in conjunction with the equilibrium curve fits 
c t 

of Equation (32). Corrector values are evaluated with averaged values of the 
coefficients after predictor values have been obtained for the numerical domain 
under consideration. 


The same general logic is applicable for the calculation of all solid 
boundaries. In particular, for the calculation of sidewalls, a local rotated 
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reference plane system is established as indicated in Figures (2) - (4). The 
calculation is performed in the plane x^ = x.^ ^ or which geometric details 
are depicted in Figure (17) for a line source or cartesian system and In 
Figure (18) for a cylindrical system. Appropriate relations In the rotated 
system are provided via the transformations 


x t " *1 
*2 * ' x 3 


(47a) 


u « u 

V a ** VJ 
W = V 


(47b) 


The metric coefficients then become 


h^ * 1 * 

h 2 « 1 

h^ » - J 2 Xj) *+ (1 - J.j + J 2 ) 


Characteristic directions in the rotated system take the form: 


h 3 dx 3 _ Hw ± a lfc 2 4w 2 -a 2 
h 1 dx 1 u 2 - a 2 


(48) 


The compatibility relation along the characteristics may be written: 
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£H_ 

BP 


d (x w/u) ± dtnP = 


-+ 
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dx 


(49) 


where 
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D ■ “ ■ [&~ u - w) A - A“ B + CJ 

h 1 h 2 h 3 ' 

* 

a - th,B ftr p- + p v- ) + j,r 2 pS - jjhjpj] 

« i * 

« 

« 

B * E (fi puv)-^ - u(h jh^pv) - - J^pw 2 ] ■ 

C » [(h^pvw)- - wO^fy^)- + J t pwuh 2 - J pmrh ] 

2 2 

The combined momentum equation takes the form: 


d© - ©f> d© -i ' dbp) 

u q u u pq 




- (Bu + Cw)J dx cl 

q SL 


where. 



_ _ [(h.h (P + pv ))- - vfh.h.pv) + J,h, (P + pw 2 )] 

hhh 1 * *9 » 3 - 2 1 

n l 2 3 ^ x 2 


The statement of the boundary condition V*n = 0 on the sidewall g = g(Xj, 
33 g(x 1 x 2 ) may be written: 



?- h 3 (g x 
u J X 1 



The system of Equations (49) - (51) is solved directly for the variables 
w/u, v/u and P. The analogous boundary relations in rotated coordinates For 
CHAR3D are detailed in Reference (4) for the three coordinate systems in- 
corporated in this analysis. 


D. Corner Point Calculations Procedure - interior corners occur in the 
internal modules and are discretely analyzed in both numerical codes by a 
redundant "weighted" characteristic procedure. For a cartesian system the 
corner results from the intersection of the surfaces z * f(x,y) and y = g(x,z) 
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Referring to Figure (19), the boundary condition V*n * 0 applied to both 
intersecting surfaces at C (the point to be calculated) yields explicit 
expressions for the flow deflection angles in the reference plane y * y c : 



tan* c 





( 52 ) 



g + f g 
tamp c = (cos^ - y 


(53) 


Then, a redundant procedure is employed wherein reference plane calculations 
for the pressure at C are performed in the reference planes z = and y = y^. 
This yields two values of pressure and Pq 7 which differ due to evaluating 
the cross derivative forcing function terms in the compatibility relations via 
backward differences. A weighting of these pressures is performed by account- 
ing for the relative wave strengths in each of these reference planes. This 
gives the stronger weighting to the calculation performed in the reference 
plane containing the dominant waves via the relation 


P, 


A %C 


+ 


(54) 


In the line source system, the upper or lower walls are specified by re- 
lations of the form z = f(r,0) while the sidewall by y = g(x,z), as indicated 

A 

in Figure (.20). Application of the boundary condition V*n at both walls in 
conjunction with the transformation 



from line source velocity components, u,v,w to cartesian components u,v,w 
yields 

,, f r + 9* tan0) + T* (g x “ tans) 
tan* = (£) = - 

C 0 + s x tan0 - g z — ) 


( 56 ) 





In’ the reference plane 0 * 0 . and the relation 

c 



g + [cos© f - sin© “] g 
x r r J z 

f A 

1 - [sin0 f + cos0 g 
* r r z 


( 57 ) 


with respect to a cartesian system (y = y ). 

c 


For a cylindrical system, the upper or lower wall is specified by aru-equa- 
tion of the form r ° f(x,0) while the sidewall by 0 = g(x,r). Expression of the 
boundary conditon V*n = 0 on both surfaces at point C results in 


and 



f 
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1 
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fe 
r r 


( 58 ) 



+ f g 
x r 



( 59 ) 


In all systems, the pressure P is determined via the weighted character- 
istic technique, while streamline relations, applied along the corner stream- 
line CD, are employed to generate remaining flow variables. 


E„ Calculation of Three-Dimensional Surfaces of Discontinuity - The nu- 
merical models described treat the plume interface (separating the nozzle ex- 
haust flow from the vehicle external flow) and plume generated external shock 
wave as discrete surfaces of discontinuity. The overall approach closely fol- 
lows that of Reference (4) and thus is summarized briefly in this report. 

The Orientation of a surface of discontinuity with respect to the x^, x^, 

Xj coordinate system is depicted in Figure (21), The local surface orientation 
is specified by an ortho-normal triad of vectors consisting of the normal to the 
surface and two surface tangent vectors. With respect to the reference plane 

A 

x 0 * constant, a tangent direction t with cosine director 6 is defined by the 

t A 

shock and reference plane intersection. The surface tangent 2. with cosine 

A A 

director a is normal to t and thus the normal to the shock n is given by the 

A A A 

cross product n = t x fc, which in terms of a and 0 may be written 


*7 



n » - eosa sing i 




sina i + cosa cosB I 

x 2 x 3 


The cross cast angle a% made Uy the cut of the discontinuity surface with the 
plane x^ = constant and the ^ coordinate direction, is related to the angles 
a and B by the relation 


tana"* = tana/cosB 


Shock Point 


For the calculation of a shock point, the local velocity vector expressed 
in the shock oriented system may be written 


- *v a, 

V = un + v £ t + 


and the Rank ine-Hugon lot jump relations in this shock normal system are: 
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Normal Momentum 
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H = constant = h+s V 
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(61) j 
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(64) j 
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( 66 ) 


(67) 


49 


•State 


P a P (P,h,$) ' (68) 

where points 1 and 2 are respectively upstream and downstream of the shack. 

Combining and rewriting Equation (63) and Equation (64) in terms of the angles 
a and 6 and velocity, expressed by q, $ and ip, one obtains 

Continuity and n Momentum 

p, u t (u 1 - \) = P 2 - P, (69) 


where 


u » q[tampsina +cosas in (3 -<Jj)] 


t Momentum 


costB-^) = q 2 cas(6-<(> 2 ) 

A , 

£ Momentum 

q^ (-s i nets i n (B — ^ ) + cosa tan^) = q 2 (~s inas in (B —( f> 2 ) + cosa tan^ 2 ) 


(70) 


(71) 


Energy 


State 


H - h + i (q 2 /cosi}; 2 ) = constant 


P 2 = P ^ P 2’ h 2’ ^ 

Referring to Figure (22), illustrating the local finite-difference network 
in the line source system, shock point C is calculated as follows: 

(l) In each reference plane J, point C is located by the relation 
h, dx, 

7— (j— = tanS. 

h i „ 

5 n 


(72) 


(73) 


(74) 


SSESSSS 
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where 3. is the shock angle in the reference plane at the preceding 
step. Having located, points x, in adjacent reference planes J - 1 
and J + 1, the cross cut angle' a" is calculated by the relation 


tana' 

c 


(x, " x« ) 

h 3 C J*H C J-1 
h^ • 2 Ax 2 


( 75 ) 


for reference planes equally spaced at intervals of Ax^. 

(2) Properties at Cj are evaluated by a standard reference plane char- 
acteristic procedure. Note that the position of C remains invariant 
In this calculation procedure and thus properties may be evaluated 
at both the predictor and corrector level. Points and (Figure 
22) are the intersections of the reference plane characteristics 
passing through with the initial data surface. 

(3) A value of the shock angle in the reference plane, 3 , is assumed, 

c 

yielding a via Equation (61). 

C " 

(4) The Rankine-Hugoniot relations (Equations 69 - 73) yield properties 
q, tji, P and p at C^. 


(5) The compatibility relation applied along the isprunning characteristic 

AgCg with $ = <J»c 2 and coefficients averaged yields an alternate value 

of pressure P . 

c 2 

(6) The angle $ c is perturbed and steps (3) and (5) are repeated until the 
pressures P C2 and ?* c are in agreement to within a specified tolerance, 
linear error extrapolations are employed to speed convergence. 


(7) Having converged in each reference plane J, the following global 
Iterative procedure may be employed although the additional ac- 
curacy provided by these additional steps has not been ascertained. 


52 


original PAGE lb 
OS POOR QUALITY 


(a) Evaluate cross derivatives S/SXg at points C ^ and incorporate 
their values into the forcing function terms of the compatibility 
relations along A 2 C 2 . Although coefficients were averaged in 
this relation, forcing function terms were evaluated based on 
cross derivative values at A 2 * 

4 

(b) Relocate the points x, replacing tanS, of Equation (74) with 

1 5 c ' 

't [tanS* + tanB ] , in each reference plane J. Then, the cross 
Z I c 

cut angles a' must also be reevaluated via Equation (75). 

(c) Repeat steps ( 2 ) through (6) in each reference plane J with 
the initial estimate of & c being the converged value from the 
first global iteration. 

Contact Point 

The contact point calculational procedure is substantially more complex 
than its two-dimensional counterpart since the streamlines passing through a 
point on either side of the contact discontinuity not only differ in the values 
of composition and stagnation properties, but also may be highly skewed with 
respect to each other. Thus, as for a solid boundary, the angle made by the 
cut of the contact surface with the reference plane differs from the stream- 
line projection onto the reference plane. Referring to Figure (23) illustrating 
the local surface geometry, the streamline passing through point C on the upper 

JU 

side (side 2) of the contact surface is the line C^ while that through the 
lower side (side 1) is the line D" C^. The angle 8 is the angle made by the 
contact surface intersection with the reference plane = constant while a' 
is the cross cut angle. Geometric relations are those presented above. A 
local iterative procedure in each reference plane, analogous to that for the 
shock point calculation is performed to satisfy the boundary conditions: 

(V»n) = (V*n) = 0 (76a) 

C 1 c 2 

P « P (76b) 

C 1 c 2 
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The computational procedure may be summarized as follows: 

(1) Points Xi are located in each reference plane J via Equation (7*0 

J c 

while cross cut angles a'' are obtained via Equation (75)- 

(2) A value of B is assumed yielding a via Equation (6l) and an ex- 
pression of the boundary condition V*n = 0 (Equation 76a) in the 
form 

sin(B - <jt ) + tana tamf; = 0 

1,2 C C 1,2 

(3) The standard solid surface calculational procedure (Section INC) 

in conjunction with Equation (77) yields values of P, <{> and \ p at points 

c t and c 0 for the assumed value of B . 

1 2 H c 

(4) The angle B c is perturbed and steps (2) and (3) are repeated until 

the pressures P c ^ and agree to within a specified tolerance. 

Linear error extrapolations are again employed to speed convergence. 

(5) Upon convergence in each reference plane J, the global iterative 
procedure employed for the shock may be employed with cross deriva- 
tives 3/3x2 now eva * ua1:ec * along grid points c^ and 

F. Cowl Lip Exhaust/External Flow Interaction - At the cowl lip, the in- 
viscid interaction between the nozzle exhaust flow and external stream is dis- 
cretely analyzed establishing the initial geometry of the contact surface 
separating these streams. For underexpanded flows, this interaction results 
in an expansion fan propagating towards the vehicle undersurface and a plume 
generated bow shock as schematized in Figure (24). The calculation of this 
interaction is simplified by recognizing that the flow phenomena is locally 
two-dimensional in planes normal to the cowl trailing edge. Figure (25) de- 
picts the coordinate system associated with this local normal plane. By 
transforming data to this coordinate system, standard two-dimensional procedures 
are employed in determining the contact angle 8 which yields a- pressure balance 








between the exhaust flow and external stream. 

* 

For a non-swept cowl trailing edge, the local unit vectors are obtained 
via the transformation 




sinB 


-cosSsina 


A = 


cosa 



cos0 -sinBsina 


trailing edge and n is normal to the discontinuity surface at the trailing 
edge. The angle a is the cross cut angle made by the cowl trailing edge with 
the Xj coordinate direction, while 3 is an angle made by the cut of the dis- 
continuity surface with the reference plane x 2 - constant. The angle 0 is the 
cut of this discontinuity surface with the local normal plane and is related 
to a and 0 by the expression 


(78) 


tan0 = cosatanB 


( 79 ) 


The iterative process is initiated by transforming the velocity components 

A 

to cowl oriented coordinates where n is identified as the unit normal to the 
inner cowl surface at the trailing edge.* The initial velocity of the exhaust 
flow in this system is expressed by 


- Tj 

V = un + vfc + wt 


( 80 ) 


where the components u, v and w are obtained via the transformation 


; For a cowl surface specified by the relation z = f(x,y), 3 is given by 

*2 * 

-1 


0 => tan (- 


r) 


1 + f 
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The Mach number projection onto the local normal plane, used to Initiate 
the Prandtl -Meyer integration, is expressed by 

,.2 ,^2 . 2 

- (u + w )/a 

Then, the interface angle B is incremented in small steps of A3, where for each step 
l the following procedure is followed: 


(1) The pressure on the interface is determined via the relation 



where 3. =6. + (i-1) A3 and subsequent projected Mach numbers are 
obtained by a standard isentropic relation. 


(81) 


(2) The pressure jump across the bow shock wave is determined in cor- 
respondence with the change in flow deflection angle in the local 
normal plane as expressed by 

6 , = cos ^ (n • n ) (82) 

sh w c 


where n is the unit normal to the outer cowl surface while n is the 
w c 

.normal to the contact surface as expressed by Equation (78) with 
•3 = 3 r 

(3) The shock angle 3 g is determined via the relation (for a perfect gas) 
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cot8 . 
sh 


tanB s 


<rM) 

2(M 2 sln 2 B -1) 
e s 


-1 


(83) 






where ?( 2 is the external stream Mach number projected onto the local 
normal reference plane as given by 

ft 2 ^ + ^ 2 )/a 2 fSIs) 

r\s f\, 

where the velocity components and are obtained via the transfor- 
mation of Equation ( 80 ) identifying 0 in matrix A with expressions of 
the normal to the outer cowl surface. 


(4) The pressure on the external side of the interface associated with 

the change in flow deflection angle 6 . and shock angle 0 g is given by 


s. 

i 


2yMg sin 2 0 s - (y-1) 

(y+D 


(85) 


(5) The process is repeated until P ^P.; then properties are obtained 

S • 1 I 

by linearly interpolating between values at the last two iterations. 


(6) Resultant velocities in the local normal system are transformed back 
to standard reference plane components via the inverse transformation 



G. External Corner - End Module Calculational Procedures - For an end 

module, the local exhaust/external flow interaction processes at the trailing 

* 

edge occur in mutally perpendicular 'planes. To best accommodate the plume 
flowfield- calculations in this vicinity, a combination of reference plane 
systems is employed. For the rectangular end module schematized in Figure (27), 
vertical reference planes (y=constant) are employed in the central region while 
horizontal reference planes (z=constant) are used in the vicinity of the module 
sidewall. A cylindrical "wraparound" reference plane system is employed in the 
region of the corner to provide for the transition between these two systems. 
Implementation of this hybrid grid system provides reference plane alignment 
essentially perpendicular to predominant wave and discontinuity surfaces as well 
as to the vehicle undersurface. 
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the details of the interaction process in the corner region ere depicted 
tn Figure (28). For the simple corner comprised of the intersection of the 
surfaces z - constant (cowl) and y = constant (sidewall) with uniform internal 
and external flowfields, the following two-dimensional regions may be identi- 
fied, for an underexpanded exhaust: 

t 

(t) A cowl interaction region resulting in an interface deflection angle 
6,j t a bow shock angle and uniform pressure (cowl) extending 
from the bow (external) shock to the terminating ray of the Prandtl- 
Meyer expansion fan emanating from the cowl. 

(2) A sidewall interaction region having corresponding interface an 

shock angles S 2 and respectively and region of uniform pressure 

V su) ' 

Both these two-dimensinal regions occur outside the domain encompassed 
by the Mach cone enamating .from the juncture of the cowl and sidewalls at 
the trailing edge. The region outside the intersection of the two Prandtl- 
Meyer fans emanating from the cowl and sidewalls remains undisturbed while the 
region within the intersection of these expansion fans is highly three-dimensional, 
numerical solutions for internal and external corners have been performed ex- 
ploiting the conical invariance of such flowfields (Reference 16). A compari- 
son of such internal corner solutions with those obtained by BIGMAC and CHAR3D, 
employing standard boundary calculations! procedures, have been demonstrated in 
Reference (12) and are described in more detail in the next section. Similar 
comparisons for the substantially more complex underexpanded exhaust/external 
flow Interaction problem pend the availability of solutions by techniques ex- 
ploiting the ocnical invariance of this flowfteld. 

Details of the localized grid network in the corner region, for a quiescent 
external stream, are provided in Figure (29). Each of the three reference 
plane systems are discretely analyzed employing overlap planes at the boundaries 
of each system to provide a mechanism for the evaluation of cross derivatives. 

The grid point located at the origin of the cylindrical wraparound network is 
calculated in both the vertical and horizontal reference plane systems and the 
results, averaged. Evaluation of the allowable marching step via the CFL sta- 
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FIGURE 29. GRID NETWORK FOR END MODULE EXHAUST PLUME IN QUIESCENT STREAM 



blllty criterion does not consider transversal spacing in the cylindrical 
system, since this would be overly restrictive in the proximity of the origin. 
Rather, the numerical domain utilized for the evaluation of such cross deriva- 
tives is "extended" to encompass grid points outside the intersection of the 
Hach four cone from the grid point being calculated with the initial data 
plane, thus ensuring stability. 

Rather limited experience in performing such corner calculations has 
failed to lead to an optimized approach for grid orientation in this region. 
However, the location of the origin of the cylindrical "wraparound" reference 
plane system cannot be fixed apriori in a generalized manner catered to a 
wide class of exhaust/external flow conditions. Rather, the origin must be a 
"floating" one, which adjusts to the local growth rate of the plume interface 
both vertically and horizontally such as to best approximate the local radius 
of curvature in the "wraparound" region. The logic to perform the above has 
been developed for Program BIGMAC, which of course includes a complete revision 
of the local corner grid network everytime the origin is revised. 

H. Multiple Module Interactions - The numerical models developed analyze 
the plume interactions associated wi th exhausts emanating from multiple nozzles 
separated by common walls, as illustrated in Figure (1). The internal (nozzle) 
flowfield calculations are performed for each individual module and the re- 
sultant exit plane data stored on local files are combined to generate a complete 
exit plane map of the exhaust flowfield. 


The resultant exit plane properties may differ in pressure and stagnation 
properties from module to module and, a discontinuity in cross flow angle may 
exist, produced by the finite trailing edge wedge angle in the common walls 
separating the individual modules. Thus, in addition to the primary interac- 
tion between the exhaust flow and external stream (as depicted in Figure 2b) a 
module to module interaction process occurs at the trailing edge of the walls 
separating adjacent modules, resulting in a wave system propagating predominantly 
in the spanwise direction. This process is schematized in Figure (30) wherein 
the wedge induced cross flow shocks interaction with the downrunning Prandtl- 
Meyer expansion fan and plume induced bow shock in the domain included within 
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the Mach cone emanating from the juncture of the wedge and cowl trailing 
edges . 

Preliminary results from calculations made in these regions indicate the 
requirement for maintaining proper grid control and calculating these local- 
ized flowfields initially employing a highly refined grid network. In parti- 
cular, the calculation of these complex interaction regions requires a finite 
number of integration steps to converge to the conically invariant solution 
whose boundary conditions are those associated with the flowfield in the immed- 
iate vicinity of the wedge/cowl trailing edge juncture. In corporating the 
calculation of such regions into the numerical codes on a grid scale commensur- 
ate with the overall flow domain may result in numerical difficulties attributed 
to the influx of waves (resulting from nonuniform flowfield properties at the 
module exit planes) into these regions before convergence has been achieved. 
Since the convergence process should take place on a length scale which is 
small in comparison to global inviscid length scales (i.e., module dimensions), 
the calculation of such interaction regimes seems most readily performed via 
the performance of a separate, localized calculation. The locally converged 
results can then be incorporated into the overall grid network as initial 
conditions made in the same manner that a Prandtl -Meyer expansion solution 
would be incorporated into a two-dimensional grid network. 


IV. SAMPLE CALCULATIONS 


A systematic procedure has been followed in assessing the capabilities 
and limitations of programs BIGMAC and CHAR3D. First, the ability of these 

I 

codes to capture shock waves and trace their propagation in the reference 
planes has been evaluated yia the performance of a series of two-dimensional 
calculations in convergent and divergent ducts. Three-dimensional capabilities 
were then evaluated via the performance of a series of corner calculations for 
which experimental results and/or conically invariant solutions were available 
for comparison. The favorable results obtained in these calculations supported 
useage of these codes in generalized three-dimensional situations. 

Results have been obta i ned with BIGMAC for the flowfield in a rectangular 
nozzle, demonstrating its ability to calculate the interactions associated with 
shocks emanating from mutually perpendicular surfaces. Applications to multiple 
module exhaust flows have indicated the requirement for 1) a refined network 
in the vicinity of the intersection of the intermodule walls and cowl trailing 
edge and 2) a floating origin for the wraparound region associated with the end 
module exhaust flow. A descripton of the calculations performed is provided 
below. 


A. Single Wedge Inlet - Calculations were performed for a 10° wedge inlet 
having the following uniform entrance conditions: 

M t = 2.94 
P 1 = 845.5 
T 1 = 2328°K 

Y = 1.4= const 

(All pressures are nondimensionalized with respect to 47.88 N/m2.) 

Results obtained for the upper and lower wall pressure distributions with CHAR3D 
are displayed in Figure (31) for 11 and 21 ooint grid networks. In performing 
these calculations, the grid point on the wedge surface was initialized by con- 
ditions behind the wedge shock. Hence, the results were shifted slightly to 
reflect the uncertainty in shock location between the first and second grid 
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points. The agrees ent with the exact solution is quite satisfactory. Shock 
locations are accurately predicted with a spread over 3 to 5 axial stations 
at the reflection points and no overshoots are encountered at the shock waves. 
The results obtained with the 21 point grid provide a somewhat sharper defini- 
tion of the shock jumps although the 11 point grid is felt to be adequate. 

Results obtained with BiGMAC utilizing a 21 point grid are depicted in 
Figure (32). The three sets of results depicted were obtained employing the 
predictor/connector algorithm in variant modes. With IFLIP^I, predictive de- 
rivatives are made with forward differences while connector derivatives are 
made with backward differences; with I FLIP- 0 the opposite procedure is fol- 
lowed; while with the FLIP-FLOP procedure the sequence is alternated at suc- 
cessive Intergration steps (i.e., IF.IP=1 at odd steps and IFLIP-0 at even 
steps) . 

Clear benefits accrue from implementation of the FLIP-FLOP procedure in 
terms of accurately locating the shock reflection points. The overall 
accuracy (with the FLIP-FLOP option) is quite comparable to that of CHAR3D. 
Shock waves appear more 'sharply defined by BIGMAC although overshoots are en- 
countered at the reflection points. These overshoots are of minimal axial 
duration and should have a negligible effect on forces and moments. 

Subsequent studies, outside of the scope of this effort (i.e, Reference 
17» Supersonic Compressor Studies) have indicated that such overshoots are 
partially attributable to the convergence of streamlines in the reference 
planes, in regions of large compressions. This convergence also leads to a 
rather restrictive forward marching step with grid points outside this com- 
pressed region being advanced with local Courant numbers substantially less 
than one. This is reflected in numerically diffusive effects whose buildup 
can become rather substantial in ducted convergent flows. Elimination of this 
behavior has been affected by incorporating grid controls in the reference 
planes; in particular, it was found tiat grid points In the reference planes 
should be dropped when they restrict the allowable marching step in comparison 
to the average allowable step by more than about sixty percent. Implementa- 
tion of this criterion has substantially improved the performance of these 





codes although the results presented herein were performed without this modifica- 
tion. . . . 


B. Double Wedge Inlet - Calculations with CHAR3D were performed for the 
double wedge inlet depicted in Figure (33), having the same uniform entrance 
conditions as the single wedge inlet case just described. Pressure profiles at 

x = .64, 1.16 and 2.14 are depicted In Figures (34), (35) and (36), respectively, 
for both 21 and 41 point grids. The entropy distribution at x = .64 with a 21 
point grid is depicted in Figure (37)- These calculations further demonstrate 
the predictive capabilities of CHAR3D for a complex two-dimensional situation 
involving multiple wall and shock interactions. 

C. Two-Dimensional Convergent Duct - Calculations were performed for the 
duct geometry depci ted in Figure (38) and tabulated below, for the same entrance 
conditions as in the previous calculations. 

' Lower Surface : 

0 < x < .5; • ? = 

-5 < x < 1 .5; z = 

1 .5 < x < 2.0; z = 

2.0 < x < 3-0; z = 

3.0 < x < 2; z « 

Upper Surface : 
x < 0, z = 1 . 0 

Results obtained with program BIGMAC employing an 11 point grid are depicted in 
Figures (38) to (41). Results are compared with those of program SEAGULL employ 
ing a 21 point grid and thought to accurately represent the exact solution. 

The calculated shock propagation pattern is depicted In Figure (38) for 
five wall reflections. The upper and lower wail pressure distributions are 
shown In Figures (39) and (4d) respectively, while the upper and lower wall 


0 

• 1(x-.5) 2 
. 1+.2(x-l . 5) 

. 2+. 2 (x-2) - . 1 (x-2. ) 2 
.3 


*Discrete shock code developed by M. Salas at NASA Langley for two-dimensional 
and axisymmetric internal flows (Reference 18). 
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FIGURE 33. DOUBLE WEDGE INLET CONFIGURATIONS 
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FIGURE 38. CONVERGENT DUCT GEOMETRY AND SHOCK PROPAGATION PATTERN 
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FIGURE 39. CONVERGENT DUCT, UPPER WALL PRESSURE D I STR I BUT I CM— B 1 GMAC 
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entropy variations are depicted in Figures 0*1 Ai and (4lB), These results 

demonstrate ability of program BIGHAC to numerically "capture 1 * a shock formed 

via the enveloping of compression waves and carry it with minimal diffusion * 

over multiple wall reflections. The favorable comparisons obtained with this 

relatively coarse 11 point grid are quite promising. "•] 

Upper and lower wall pressure distributions as obtained with CHAR3D employ- 
ing an 11 point grid are depicted in Figures (42) and (43). These comparison ■ > 

indicate that BiBMAC and CHAR3D provide results of comparable accuracy for the I 

same grid definition. It is felt that the diffusive behavior exhibited after ') 

4 or 5 wall reflections would be largely eliminated by the grid control tech- J 

nlque previously described. '] 

V \ 

i 

i 

0. Two-Dimensional Divergent Duct - Calculations with BIGHAC were per- j 

formed for the nozzle geometry depicted in Figure (44). The results are again j 

compared with those of SEAGULL. The solid lines in Figure (44) depict the j 

polynomial approximation of the nozzle surfaces utilized in the SEAGULL calcula- j 

tion, while the discrete points are those obtained with the spline fit geo- j 

metry package described in Appendix II and employed In the BIGHAC calculation. j 

It is noted that a rather poor fit for the upper watl contour was obtained in i 

the vicinity of x = 7 with the spline fit approximation, which might have been . 'j 

'1 

avoided by utilizing more contour data points in generating the spline fits in j 

this region of rapidly changing curvature (see Reference 11)* i 

j 

j 

Initial conditions for this calculation were again the same as in previous * 

cases* Resultant wall pressure distributions obtained with BIGHAC are compared 1 

with those of SEAGULL for the lower and upper nozzle walls in Figures (45) and j 

(46) respectively* The poor geometric fit for the upper wall in the vicinity j 

of x = 7 is clearly reflected in the oscillation depicted in Figure (46) in .j 

this area. The deviations in upper wall pressures downstream of x - 12 are j 

not readily accounted for although the spurious waves emanating from the upper ] 

wall around x = 7 may have a dispersing effect on the expansion waves emana- ■; ] 

ting from the lower surf ace, contr i but i ng to this beha ior. Radial pressure 
profiles are compared at the axial locations x ■v 2*3, x 5.1, x ^ 9 and 
x * 13-5 in Figures (47A,B t C and D) respectively, indicating favorable agreement. 
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FIGURE 47. DIVERGENT DUCT, PRESSURE PROFILES 






E. internal Corner Calcualtions - Comer flowflelds associated with 
interacting waves formed from, mutually perpendicular surfaces represent ideal 
cases with which to test the developed codes. Experimental data and calcula- 
tions exploiting the conical invariance of the reported cases is available for 

ft * 

comparison. In all cases reported, flow at the initial station was uniform 
and the interaction flowfield was generated by an abrupt change inwall angle 
at the initial station for two mutually perpendicular surfaces. 

Results for a 5° double expansion corner are depicted in Figures (48) and 

(49) . These results were obtained with CHAR3D starting from the conditions 
P * 945.5 and M =* 2.94 for an 11 x 11 grid in cartesian coordinates. Re- 
suits are depicted after nine axial marching steps as required to establish a 
converged solution. This is indicated by the axial variation or corner pressure 
which originally overexpands and then recompresses to the converged value. 

The finite convergence length required to achieve invariance has important 
Implications in the generalized application of corner boundary procedures in 
regions of discontinuity. In particular, one must achieve convergence in a 
distance which is small compared to the overall inviscid length scale of the 
problem - this sets the grid network size for the localized corner calculation, 
tn addition, one must isolate the computation of this region so that extraneous 
waves doe not interfere with the convergence process. 

Results for an expansion-compression corner generated by an abrupt 5° ex- 
pansion and 74° compression of a uniform Mach 2 flow are depicted in Figure 

(50) as calculated by B1GMAC. An 11 x Vi cartesian grid was employed and the 
converged results are depicted after 10 axial marching steps. These results 
are compared with the detailed conically invariant numerical solution of Shankar 
(Reference 19) and the experimental results of Nangia (Reference 20). 

Results for a 12.2° double compression corner are presented in Figure (51) 
as calculated with BIGMAC. A line source reference plane network was employed 
in this calculation in view of the larger turning angles than in the previous 
two corner flows. The initial flow was at a Mach number of 3*17. Converged 
results are depicted after 35 axial marching steps indicating an increase in 
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•convergence length with the severity of the discontinuity as would generally 
be expected. Comparisons are made with the numerical solution of Shankar 
(Reference 19) and the experimental results of Charwat and Redekeopp (Refer- 
ence 21). 


F. Square Nozzle - The three-dimensional flowfield within the square 
nozzle depicted in Figure (52) has been calculated by B1GMAC. This flowfield 
is characterized by the initial interactions of expansion waves emanating from 
mutually perpendicular surfaces cud the subsequent interaction of enveloping 
shock systems generated by reccmpression on the upper wall and sidewall. The 
initial Mach number was 2. 94 and the initial nondimensional pressure was 845.5. A 
perfect gas calculation with y = 1.4 was performed. The calculation employed 

21 grid points in each reference plane with 11 reference planes initially 
(reference plane number 1 was the plane of symmetry). A cartesian network 
was utilized and additional reference planes were inserted as the module side- 
wall opened. At the straight section, the network contained 18 reference 
planes. Pressure contours on the symmetry plane y = 0 are depicted in Figure 
(53). It is of interest to note that the contours on the symmetry plane z = 0 
are virtually identical to those of Figure (53) thus providing a check on the 
overall symmetry of the computational system. Similar checks with CHAR3D have 
not provided the required symmetry in cases where strong wave systems were 
propagating normal to the reference planes. Of particular interest in Figure 
(53) is the intersection of four three-dimensional shock surfaces at x * 17 and 
and y “ z * 0. This results from the reflection of the envelope shock produced 
by the sidewall and the reflection of the envelope shock produced by the upper 
wall, resulting in an approximate 15/1 pressure ratio at this location. The 
axial pressure variation along the corner is depicted in Figure (54) while pres- 
sure variations along several streamlines in the symmetry plane are depicted in 
Figure (55). 

G. Single Module Nozzle-Exhaust Flowfield - The internal nozzle and ex- 
haust plume flowfield associated with the single module depicted in Figure (56) 
has been calculated by BIGMAC. The gas mixture was assumed to be perfect and 
the calculation was initiated by the following urt'form conditions at the com- 
bustor exit: 


97 








26 k l! GRID 'Mgfg.94 



Figure 5**. Streamline pressure distribution at sidewall corner of square nozzle. 
* -21 *11 GRID M^-2.94 



Figure 55* Streamline pressure distribution in plane of symmetry of square nozzle. 
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TS*e external flow was quiescent with, v = 1.4. 

* 00 


The Internal flowfleld was two-dimensional and calculated employing an It 
point grid up to the cowl exit station (x ■ 2.98). Note that x is measured 
along the vehicle undersurface and y normal to it. The resultant flowfield at 
the cowl exit station is underexpanded along the cowl trailing edge and mixed 
(partially overexpanded and partially underexpanded) along the end wall trail- 
ing edge. The reference plane networks employed are sketched in Figure (57} - 
Fifteen vertical reference planes, fifteen cylindrical reference planes in the 
corner wraparound system and four horizontal reference planes are employed in 
this calculation. Resultant interface locations are depicted in Figure (58). 
The last marching station calculated was at x/H t * 6.43. The requirement for 
relocating the cylindrical origin (floating origin option) is clearly necessary 
to continue the calculations beyond this point. All cases with some degree of 
overexpansion will require this modi ficiat ion to account for the nonuniform 
collapse .in plane size. The growth of the plume Interface in the plane y/H t *3 
Is depicted in Figure (59)* while Isobars at the stations x/H fe *3»04, 3*75* ^-99 
and 6.43 are depicted in Figures (60) - (63) respectively. 

H. Double Module Nozzle-Exhaust Flowfield - For this case, the side 
view is the same as that depicted in Figure (56)» while the top view is illu- 
strated In Figure (64). Initial conditions at the combustor exit are the same 


1Q2 


ORIGINAL PAGE lb 
DE POOR QUALITY 










: mm r 

L 1 



FIGURE 




as for the single module case except for the pressure ratio P^/P^ 

which is now 10.5 . Resultant plume shapes are de- 

plctedf in Figure (65). While the ability to calculate such complex flowfJelds 
Is deiRsnstrated, the results indicate the requirement for major modification* 
In grid control procedures. 


* It should be noted that this flowfield was calculated on a one-shot basis. 
Program B1GMAC calculated the internal flowfields for each of these modules 
employing a line source reference plane system up to the cowl trailing edge 
station (x/2 » 3.7) and automatically Interpolated the exit plane results 

b 

for the two modules in a Cartesian framework and proceeded with the calcula- 
tion of the exhaust plume flowfield. The initial angle of the contact sur- 
face separating the exhaust flow from the quiescent stream Is explicitly 
calculated as described in Section I IIP. The intermediate interaction pro- 
cess is also explicitly calculated. The locally two-dimensional Interaction 
pressure and flow deflection are computed on either side of the module 
juncture. En'addition, at the module juncture-cowl Intersection solution 
computed above. Thus, in this test case the cowl plume at the module junc- 
ture experiences a diminished level of overexpansion and possibly a slight 
underexpansion as a result of the module Interaction. At present, in lieu 
of an exact conical solution, at the cowl module juncture, the procedure 
described above appears acceptable. However, detailed asymptotical conical 
solutions for this problem should be obtained. 

An interesting comparison is made in Figure (66) employing a discrete 
floating shock technique as reported in Reference (22). The calculations are 
for the impingement region of two Initially uniform rectangular jets ex- 
panding to background conditions. The comparison, which Is meant to be 
qualitative , shows strikingly similar plume shapes, including the plume 
kink and sharp peak. The results of Reference (22), while preliminary, do 
lend credibility to the intermodule interaction flowfield as computed herein. 


Pressure contours at the cowl trailing edge station {x/Z^ * 3.7) as well 
as the axial stations x/Z fc - 4.77, 6.1A and:‘9.93 are presented in Figures 
(67) " (.70)., respectively. 
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V. CONCLUDING REMARKS 


Programs BIGMAC and CHAR3D have been developed to provide predictive 
computational tools for the analysis of supersonic three-dimensional nozzle* 
exhaust flowflelds. Preliminary applications of these codes to a variety of 
two and three-dimensional situations indicates that the performance level Is 
quite satisfactory for this purpose. Both codes have demonstrated the ability 
to capture shock waves in the reference plane and to predict the shocks propa- 
gation pattern even after multiple wall and shock Interactions. Program BIGMAC 
has demonstrated this ability in complex three-dimensional situations wherein 
the Interaction of shock waves from mutually perpendicular boundaries has been 
calculated. 

A comparison of these two codes indicates that BIGMAC provides non- 
preferential treatment of general three-dimensional flows while CHAR3D is 
better catered to flows with wave systems propagating predominantly in the ref- 
erence planes." Difficulties have been encountered with CHAR3D in attempts to 
analyze flows with strong shock systems travelling normal to the reference 
planes. This is attributable to the basic reference plane characteristic ap- 
proach wherein cross flow derivatives are treated as forcing functions in the 
Integration step. If these forcing function terms are large and changing 
rapidly along the streamlines (as is the case for cross flow shocks), the 
approach tends to breakdown. For such flows, BIGMAC is clearly the prefer- 
able code 

The effectiveness of these codes is largely attributed to the treatment 
of the boundaries. All boundary points are analyzed by a characteristic pro- 
cedure in reference planes which attempts to include the local boundary normal. 
For corner flows, a redundant approach employing two perpendicular sets of 
reference planes has proven quite effective. 

Initial attempts at analyzing the flowfields downstream of the cowl exit 
plane has led to the requirement for several modifications to the current 
approach. In particular, complex conical interactions occur at the juncture 
of mutually perpendicular trailing edges. It is suggested that these localized 
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Interaction flowflelds be Individually -solved by subroutines bmilt fwto 
these codes* The converged solutions way then b* Incorporated into tfe* &mr- 
all ststerlcal grid system, thus eliminating the requirement for an overly 
refined grid, network In these Interaction regions. In addition, a ^floating” 
cylindrical origin appears necessary for the wraparound network tea the vicin- 
ity of the end module exhaust, and additional grid control techniques are 
required In the region downstream of the Intermodule ^seckje trailing adfces. 

The success obtained In the preliminary calculations par formed to d©te 
clearly supports useage of these codes in future studies. The models developed 
are by no means limited to the calculation of nozzle-exhaust flows. Program 
SUPFA& (Reference 17) is a revised edition of Program Bi&tftC which calculates 
the three-dimensional flowfield equations in a supersonic fan stage. Similar 
modifications can extend the applicability of these models to more generalized 
three-dimensional situations. 
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APPENDIX A 


CURVE FITS FOR T, h and p 

The variation of r (the equilibrium value of y) as a function of tempers" 
ture (T), pressure (P) and equivalence ratio (4) is presented graphically In 
Figures (At-, A2, & A3) , from values tabulated in Reference (4) . in Figure (At) 
It can be seen that F is a strong function of T over the temperature range of 
Interest, while the effect of varying composition is small by comparison. More- 
over, Figure (A'T) indicates that r is moderately sensitive to pressure and the 
degree of sensitivity Increases substantially as the temperature level increases 
and dissociation effects become important. 

As a result of these observations, temperature Is the primary independent 
variable, while pressure is the secondary independent variable and composition 
acts as a perturbation variable. Thus, we can fit the function r (T,P,4) with a 
polynomial in T. and add on a temperature dependent correction term for the ef- 
fect of pressure and a temperature independent correction terni for the effect 
of 4. * - 

An examination of Figure (Al) suggest that the function r (T) can best be 
curve fit by breaking up the temperature range into three intervals such that 

5 

the function can be represented by a parabola In each range. Choosing p “ 10 
pascal and 4 * 1 as our base we, therefore, find three functions 


r^T.IO 5 ,!) - - 1.833xU»“ 7 T 2 + 7*5x10" 5 T + 1.3&7 
r 2 (T,10 5 ,l) * 2.0x10" 8 T 2 - 1 .38x10”** T + 1.423 
r 3 (T,1Q 5 ,1) =■ 7*27x10 " 8 T 2 - 4.57x10“** T + 1.85 
and define the basic temperature function as 



^(TJO 5 ,!) 


0 \ 

T < 500°K 

r{T,io 5 ,i) - { 

r 2 (T,io 5 ,D 

^ for < 

500 <_ T <_ 2000°K> 

1 

r 3 <T,io 5 ,0 


T > 2000°K 

V *“ 


( 1 ) 

( 2 ) 

( 3 ) 


( 4 ) 


tig 
















figure (A3): Indicates that ||- Is constant In the two ranges 4 < 1 and 4 > 1, 
but Is a function of T ? Fitting the function ^ in each of the ranges of $ 
we obtain 



for 

*« < r 


* 

4 :> 1 


where 

n^T) » 4x1 Cf 9 T 2 - 2x1 0~ 9 T - 0.019 
n 2 (T) » 3.39x10~ 2 T 0 * 5 - 3.91x10"** T - 0.681 


This now defines F as a function of both temperature and 4 by means of the 
equation 


r(T,io 5 ,4) * r(t,io 5 ,i) + (* - i) U 


( 5 ) 

( 6 ) 
(7) 


( 8 ) 


Finally, the effect of pressure must be Included. From Figure { 18) we observe 
that r may be approximated -as 


r(T,P,t) - r(T,10 5 ,*) + m [log I0 (px10 5 ) - 5] 


( 9 ) 


where m Is a function of T. Deriving m, we find 


m 


-2.15x10" 8 T 2 + 0.91x10”^ T - 0.0695 


1 


for 


Summarizing, the final function obtained is 


T < 1000 K 


• 1 T > 1000 K 

'J 


( 10 ) 


r(T,P,4) * r(T,io 5 ,l) + m(|2_E - 5 ) + §£• (4 - 1) 

where the functions r(T,10 9 ,1), and m are given by Equations (4), (5) and 
(10) respectively. 


(It) 


The curve fit for enthalpy is derived in a similar way. Figures (A4) & (A5) 
present the variation of h with temperature, pressure and equivalence ratio. As 
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FIGURE A5. enthalpy as a function of pressure 
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a quadratic function of T, 
additive terms for the ef- 
irfzed below. 


for T > 2000°K and $ 1 

a ? « 10“ 6 (1.792* 2 V.3983* + .310) 

» 10‘ 3 (-9.05$ 2 - .0791 7$ * .245) 

C| * 10. 86* 2 - .11834) + .970 

for T > 2000°K and $ > 1 

Bj - 10 -6 (4.8l$ 2 - 13.94* + 11-59) 

bf » 10” 3 (-23-08$ 2 + 66.82$ - 52.61) 

C| 27.05$ 2 - 73.73$ + 58.39 

When the Inverse function T(h,$,p) Is required, it is obtained by an iterative 
solution of Equations (12) through (18). 

The density Is found by obtaining a curve fit for the mixture molecular 
weight and using the equation of state 

p -fiS 1 

RT 

where R is the universal gas constant and m is the molecular weight. 


(17) 


if 


m 


08 ) j|| 


: A 

IM 


i' .1 

fi 




iKJ 


(19) 
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Hie behavior of m with T, p and $ is illustrated in Figures (A6) S (A7). 
We see that for temperatures less than 2000°K, m is essentially independent of 
temperature. The discontinuity in slope of m($) shown in Figure (A7) requires 
that the equivalence ratio range be split in two. Thus, 


T < 2000°K 


m($) * 



> 


1 

- 5.8954* + 28.955 

f for i 

$ <1 

/ 

- 10.6$ + 33.6 

I 

$ > 1 

i 


) 



1 1 


(20) ! 
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MOLECULAR HEIGHT AS A FUNCTION OF TEMFERATURE AND PRESSURE. 



FIGURE A7. MOLECULAR WEIGHT AS A FUNCTION OF EQUIVALENCE RATIO FOR T<2000°K 



For the higher temperature range, It is convenient to employ the form 
m ® m(*) - 6{p,$»T) 


where 


6 


C 


T-2000, 
1000 } 


n 2 4) 


and 


a, 


(4W" + b, (?*£) + 


1.5 


Mi 


2 2.3 


2 '2.3 


for 


0 < * < 1 


and , for 


1 .< * < 2 


a 2 « -2.3* 2 + 4. OH + 1.736 
& 2 » 8.61* 2 - 15.42* - 6.66 
c 2 - -16. 88* 2 + 33.21* + 14.58 
n 2 = .4375* 2 + .0625* + 2.08 

a 2 - -.822* 2 + 2.363* + 1.905 
b 2 = 2. 76* 2 - 7.56* - 8.68 
c 2 =* 3. 6* 2 + 7 . 36 * + 27.15 
n 2 - .47* 2 + 1.825* + .350 


( 21 ) 

( 22 ) 

(23) 


(24) 


(25) 


ssss ® 5 


APPENDIX B 


THREE DIMENSIONAL SURFACE REPRESENTATION AND INTERPOLATION PROCEDURES 


Consider the three-dimensional continuously differentiable surface 
z “ f(x,y) depicted in Figure (Bl), prescribed by contour data in the form 
(yjfZj)j at discrete values of x.. Prescription of the data in this form 
Is generally obtainable for most aerodynamic bodies and greatly simplifies 
the chore of numerically fitting the surface by reducing the problem to 
the determination of one-dimensional partial cubic splines in two coordinate 
dt rections. 

Assume that there are J(i) contour data pairs (y.,Zj)., (1 < j ^J(i)) 
at I contour stations x., (1 i <_ I). The number of contour pairs used to 
determine the surface and their relative spacing is arbitrary, as is the 
spacing between contour stations. We seek to determine a surface fit 
z « F(x,y) that will yield accurate values of the unit surface normal 


N 


Wv 


C1+(F X )* 


'» - (F v>* 'v 

+ (F ) 2 )* 

' y x'- 


0 ) 


In addition to prescribing the discrete contour data, and conditions 
must be specified at the bounding surface curves. A, B, C and D. AT bound- 
aries A 6 B, one would generally stipulate: 


(!y) x at -i * 1 or for ' = 1 1 2 * »l (2) 

or partial second derivative ratios 



At boundaries C & D, one would generally stipulate: 



y 


at 1 


1 or I for j * 1 , 2, — J(i) 


( 4 ) 
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MENS I ONAL SURFACE DESCRIBED 8Y 
SCRETE CONTOUR DATA 


or partial second derivative ratios 



(t Is advantageous to map the surface onto one having rectangular 
boundaries when projected onto a plane of constant z. This Is accomplished 
by the transformation 

x « x 


n * 


y ~ y A(x) 

y B(x)" y A(x) 


( 6 ) 


Z « 2 


Then, In (x,n,z) coordinates, the surface is bounded by n=o (A), n"l (B) » 
x n Xj (C) and X*^ (D) . The analysis may be extended to sweptback 
surfaces as depicted in Figure (&2) which are mapped onto rectangular boundaries 
when projected onto a plane of constant Z by the transformation 

e - !l x . c .(y-) — 

X 0(y)“ x C(yl 

(7) 

y " y A(x) 
y B(x)’ y A (x) 


2 * Z 

Such a surface fit requires specification of contour data on sweptback 
contours 5 - constant as shown in Figure (B2) (or alternately on n * constant 
if this proves more practicable as would be the case for the wing surface 
of Figure (S3)* 


okbinai- page to 

OE POOR 


133 



For simplicity, we will assume that contour data is prescribed at 
stations x « constant and we are working In the x,n>z coordinate system 
.given by Equation (6). We require an ordered array of coefficients 

(which will consist of partial second derivatives) to fit the surface 
z - F(x,y), analagous to the coefficients Mj of Section li used to fit 
the curve y(x). 

The numerical procedures entailed in obtaining these coefficients 
are as follows: 


(a) W© require spline fits for the boundary curves A S B In the 

form and using discrete data pairs (x^y^j for A and 

C*{fV(I j(|)) ^ or ® and 3 stipulations of end conditions. We obtain 

these fits employing the procedure outlined in Section II, obtaining 

coefficients M. and H B where 
A i B i 

M A| “ y A(x.) and V " y B (x j ) 

(b) We obtain fits of z(n) at all contour stations Xj for (1-1,2,— , l) 

using the contour pairs (y. . » z. .). This is done as follows for a 
■ i »J * .J 

given station x. 


(1) Obtain the contour pairs in the transformed system, 
namely (n, , , z. .) employing the transformation given 
by Equation (6). 

(2) Transform the end condi tons at A & B. if first derivatives 
were sti pul a ted, then we would require (Sz/Bn)^ which is 
simply 


Isn) * ” ; * V B{xj> "" V A(x.)^ 

1 »J 1 

(I«1 or l) 


(9 ) 
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If a ratio of second derivatives was stipulated at the end points, 
this ratio remains unchanged In the transformed system. 

(3) .We obtain the coefficients Mn- • employing the one-dimensional 

• » J 

method of Section 1 1 where 




-x. 


(32) 


(c) We now wish to obtain values of the coefficients H on a new grid spac- 
ing such that the spacing in the variable is the same for every contour 
(independent of i) as shown in Figure(B^) although not necessarily evenly 
spaced. Let n^ denote the specified n grid for k=l,2 — — -K. Then, values 
of the dependent variable z are obtained using the relation 




i,j-l 6&ru - 


h, . (l wi -' )3 

6&n. . 

I.J 


A 2 

* (z. . . - M ^-i,j) 

',J-1 n. j., 6 


n i \~ n k 

(. 1 tX Si) 


&n 


‘ »J 


* «*l.j -1 -V . , ^ 

,J i , j -1 6 An. • 

« »J 


(33) 


where n,^, < n k < n (jJ and An-.j * n. j - n 1>JM 


The coefficients M on the i,n grid are given by 


H 


n 


i.k 



i _ n. 

(_LLJL) + m 


An 


n 


I.J 


I.J 


An i,j 


(3 1 *) 
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(d) W© new obtain fits of Z(x). along the lines n^ * constant for 

£k » t,2 - — — ,K) using the contour pairs (xj* Zj j,). For n^ “ constant 

this- ts done as follows: 


Transform the, end conditions at C and D. If first derivatives 
were specified, then the derivative (z^)^ is- required under which 
j obtained employing the relation: 


(15a 

i,k 


(“) 

Wy 


( 22 .) 

v 3n'x 


i.k 


I.k 


j^ACx) j \ Y B(x) (| 3> 

Y e(x) - y a(x) 


where (|§) Is specified at I « lor I and (“) is obtained from 
£>x y dn x 

the previous spline fit z(n) at x *» const, using the relation 


(~) 

Wx 


An 


M. 


i.k 


i ,k 


An. 


- (M 


'k 


, > 

I ,k %k-1 


A \ 

TT" 




where An^ “ n^ - n k _^ 

If a ratio of second derivatives was stipulated at the end points, it 
Is assumed not to change in the transformation. 


(2} Vie obtain the coefficients M 
Section tl where 


>„ employing the methods outlined in 
*i,k 


M 


x i,k 


(^f) 

3x 


(is) 


n 

I.k 


The techniques outlined for obtaining the coefficients M n and M on an 

i,k x : j, 

ordered i ,k array are readily extended to surface fits in cylindrical coordi- 
nates where ft is desired to approximate the surface 


r » G(x,e) 


(16) 
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Assuming that contour data is given in cartesian coordinates as data pairs 
of the form (yj, 2j) ( at discrete stations x { as, depicted' In Figure (05) • 
We or course assume that the axis of the cylindrical system Is parallel to 
the x axis of the cartesian system (the axis of the cylindrical system Is 
given by y * y* # x « z*). Then, the input contour data Is converted to 
cylindrical coordinates by the transformations 

x«x 

r » ((y-y*) 2 + (z-z*) 2 )* 

0 ■ tan" 1 ((z~z*)/(y-y*) } 

Then, we perform the transformation 


x ■ x _ 
<** 

r ■ p 


( 18 ) 



139 



Referring to Figure (4), the surface z » F(x,y) is numerically repre- 
sented by the arrays z. . , M- , M x< , x. , for i - 1,2 — , 1 and 

• » K *i,k i ,k 1 K 

k « 1,2, ,K as well as the arrays Y A j, Y B ., M ft . , Mg. f° r • ** 1*2 “# 1* 

The arrays are determined by the procedure outlined in Section Mi. Hew we 
are given arrays and values of the independent variables x and y and seek to 
determine the independent variable z as well as the partial first derivatives- 
(zk) and (zy) , which suffice in determining the surface unit normal. The 

y x 

details of this procedure are as follows: 


(a) Determine Y A(x) * Y B(x)’ Y A(x) * Y B(x) the relations 


K. (V*) + 

B* ^ 6Ax. B* 6Ax. 

i i 


5w“ A - 


4 * * J* 6x. X-K. * 

* (V - \ -J-) (— 5^4 

B l ' B 1 A Ax. 


(19) 


B 1 1 


.. Ax. x.-x 

— !-) M-> 

B 1 6 ix, 


and 


Y fl M. 

Sw “ ‘ S" 


( x i~ x ) + m a r^i-D 

B ? 2 Ax. 

I 


/X~X. 


2A x. 


+ (V* Y a ) J 


A. A . . , 

B 1 B r Ax i 


< M A. - % > 

B‘ B 1 ' 1 6 


( 20 ) 
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where x M < x < x, and Ax, e x,r« r .“ 

(hj Evaluate n(x,y) employing Equation (6h 

(e) Ascertain the local grid In which the point x,n falls 

^1-] ^ ^ ^ Xj 

V? ^ 71 i \ 

as depicted In Figure <B6). 


(d) Determine z , (2 ) : 2l t ( z \ 

a nn # x ’ V 1 nn^x * 

*e* k**) « * z j» (2 ) 

c xx n d* xx'.n 


alloying the relations: 


* >N ' +H 

b *1-1. jc-> ^ ^ 




k-1 - M 
k 


** | \ (x,-x) 

i x T-)^ — 

1-1. jj-i ' 1 




b_ t * M 

k-I x 




i x i \ (x - x i-i> 




SL " 

C Hj 

i \ 


•‘V>* 


1-1, k-l 64 \ 


(*!•'• k 'V., 


1 

.3 

1 

.3 

JT 

1 

n l-U 1 

k 6i \ 
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Ab 2 
An k v 

(n k -n) 

”g“ j 

tr 

<1 

k \ («• 

Vi 1 


T -K 
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(z ) - M 

™ * Vl , k-I 

£> k 

x»const 


(x,-x) 


Ax, 


+ H 
• n 




AXj 


m 


(2 ) « M 

** C *1-1, k-l 
d I 

n**const 


V" 1 

Att,. 


+ H 


tn 'Vi> 


1 - 1 . k 4 "k 


( 25 ) 


where AXj 00 x, - Xj_^ 


“\ ■ VV 

(e) Determine z(x,y) 
z(x t y) * I ^ 


(n k - n)^ 


nn a 


+ Z 


nriL 


• Vi 1 
6 An. 


4n k (\- "> 


* (z, - z —*■£-) — — 
3 nn a 6 An. 


An \ . (n ‘ Vl } 


+ £* h - z —A~) ' 

b nn. 6 An. ! 


nn b T 


(26) 


♦ if z <x i ‘ 
2 xx — — 
c nAx . 
«- 1 


+ z 


Ax? (x. - x) 

+ (z - z — L_) ! 

c 'xx, 6 1 Ax. 


(X - Xj) 3 

xx d bAx. 


+ (z - z N (X ' X l-' ) 
U 6 Z XX J T-J -ST — 


J 


ORIGINAL PAGK & 
OF POOR QOALny 


m 


(f) Determine (2 )• 

(rv-n) 

{2 ) as — j ■ - - 

n x nn a ' "nn b 2An k 


+ z 


(n - n k _,r 


-'O 


b 3 ' /_ _ V “"k 


An. 


+ 7 . - {z " z ) —r— 

nn k nn ' ° 


Ani, 


( 2 ?) 


(g) Determine (z ) 
x n 

2 

(x. - x) (x - X. -) 

( 2 ) . - z 1 + z 11 

v x n xx 2Ax. XX. 2 Ax. 

ci a 1 

* (2 d ~ 2 c> , . , to . 

Ax. ZjtX H Zj<X "5” 1 

1 a c 


(28) 


(h) Determine (z ) 

* y x 


^ z y^x " ! ^ z n^x / ^ y B(x) " y A(xP 


(29) 


(I) Determine (z ) 
a y 


(z x>y * t2 x>n ’ ( Vx [ (, ' n)y A(x) +ny B(x)l 


(30) 
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